import matplotlib.pyplot as plt
from skimage import data, draw, color, transform, feature, measure, morphology
import SimpleITK as sitk
import cv2
import numpy as np
import time
from mytest.mycode import geoutils
import math


def circle_detect(image):
    # 加载图片，转换成灰度图并检测边缘
    # image_gray = color.rgb2gray(image)
    image_gray = image
    edges = feature.canny(image_gray, sigma=2.0)
    print(edges.shape)
    fig2, (ax1, ax2) = plt.subplots(ncols=2, nrows=1, figsize=(8, 4))

    ax1.set_title('Original picture')
    ax1.imshow(image)

    ax2.set_title('Edge (white) and result (red)')
    ax2.imshow(edges)

    plt.show()

    # TODO: 霍夫椭圆检测似乎不好用在这里？
    # # 执行椭圆变换
    # result = transform.hough_ellipse(edges, accuracy=20, threshold=250, min_size=100, max_size=120)
    # print(result)
    # result.sort(order='accumulator')  # 根据累加器排序
    #
    # # 估计椭圆参数
    # best = list(result[-1])  # 排完序后取最后一个
    # yc, xc, a, b = [int(round(x)) for x in best[1:5]]
    # orientation = best[5]
    #
    # # 在原图上画出椭圆
    # cy, cx = draw.ellipse_perimeter(yc, xc, a, b, orientation)
    # image[cy, cx] = (0, 0, 255)  # 在原图中用蓝色表示检测出的椭圆
    #
    # # 分别用白色表示canny边缘，用红色表示检测出的椭圆，进行对比
    # edges = color.gray2rgb(edges)
    # edges[cy, cx] = (250, 0, 0)
    #
    # ax2.set_title('Edge (white) and result (red)')
    # ax2.imshow(edges)
    # plt.show()


def gradient_test():
    test_image = [
        [0, 0, 0, 0, 0],
        [0, 1, 1, 0, 0],
        [0, 0, 1, 1, 0],
        [0, 0, 0, 1, 1]
    ]
    test_image = np.asarray(test_image, dtype=np.float32)
    dx_image = cv2.Scharr(test_image, cv2.CV_32F, 1, 0, borderType=cv2.BORDER_REPLICATE)
    # print(dx_image)
    dy_image = cv2.Scharr(test_image, cv2.CV_32F, 0, 1, borderType=cv2.BORDER_REPLICATE)

    # print the gradient of each vertex
    # for y in range(4):
    #     temp = []
    #     for x in range(5):
    #         temp.append((dx_image[y, x], dy_image[y, x]))
    #     print(temp)

    # calculate the tangent of [2, 2]
    dx = dx_image[2, 2]
    dy = dy_image[2, 2]
    k = -dx / dy
    print(k)

    # print("gradient of vertex lt: ", (dx_image[1, 1], dy_image[1, 1]))
    # print(np.asarray(grad_image))


def func():
    l = np.asarray([[1, 2], [3, 4]])
    m = np.asarray([[5, 6], [7, 8]])

    return (l[0], m[0]), (l[1], m[1])


def func2():
    a, b = 1, 2

    def ff():
        return a + 5

    a = 2

    return [ff(), b]


def func3(a=3):
    print(a)


def fff(array):
    array = array.swapaxes(0, 1)


def ttt(path):
    with open(path, 'w') as f:
        f.write(str(0.23123004))


if __name__ == '__main__':
    # label_path = "../after_post_process1/prediction_00177.nii.gz"
    # labelImage = sitk.ReadImage(label_path)

    # plt.imshow(sitk.GetArrayFromImage(labelImage)[140])
    # plt.show()

    # labelArray = sitk.GetArrayFromImage(labelImage)[140]

    # circle_detect(labelArray)

    # gradient_test()

    # ll, mm = func()

    # print(ll[0].tolist(), ll[1])
    # print(mm[0], mm[1])
    # lll = np.array([
    #     [[1, 2], [3, 4], [5, 6]],
    #     [[1, 2], [3, 4], [5, 6]],
    # ])
    # print(lll, lll.shape)
    # ll = lll.swapaxes(0, 1)
    # print(ll, ll.shape)
    # print(lll)

    # a = np.asarray([
    #     [[1, 2],
    #     [3, 4]],
    #     [[5, 6],
    #      [7, 8]]
    # ])
    # aa = np.flip(a, 2)
    # print(aa)

    # start = 191
    # end = 368
    # length = end - start
    # temp = start
    # print('method1')
    # start += length // 3
    # end -= length // 3
    # print(start, end)
    # print('method2')
    # start = temp + round(1/3 * length)
    # end = temp + round(2/3 * length)
    # print(start, end)

    # v0 = np.asarray([[1., 1., 0.], [0., 0., 0.]])
    # v1 = np.asarray([[2., 2., 2.], [0., 1., 0.]])
    # np.logical_or(v0, v1, out=v0)
    # print(v0)
    # print(v0.any(1))

    # a = np.asarray([[
    #     [1, 0, 0],
    #     [0, 0, 0]],
    #     [[0, 1, 0],
    #     [0, 0, 0]],
    #     [[1, 0, 0],
    #      [1, 1, 1]]
    # ])
    # print(a.any((1, 2)))

    # a = [1, 2]
    # a[0], a[1] = a[1], a[0]
    # print(a)
    # s = '(-69.5'
    # a, b = s.split('.')
    # print(a, b)

    # image_array = np.asarray([
    #     [0, 0, 0, 0, 0],
    #     [0, 0, 1, 1, 1],
    #     [0, 1, 0, 1, 0],
    #     [1, 1, 1, 0, 0],
    #     [0, 0, 0, 0, 0]
    # ])

    # array = sitk.GetImageFromArray(image_array)
    # image_array = image_array.astype(np.uint8)
    # image_array = image_array[:, :, np.newaxis]
    # contours, hierarchy = cv2.findContours(image_array, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    # print(contours[0].shape, hierarchy)
    # rect = cv2.minAreaRect(contours[0])
    # print(rect)
    # print(cv2.boxPoints(rect))

    # rad = np.deg2rad(90)
    #
    # R = np.asarray([
    #     [np.cos(rad), -np.sin(rad), 0],
    #     [np.sin(rad), np.cos(rad), 0],
    #     [0, 0, 0]
    # ])
    # v = np.asarray([1, 0, 0])
    # # v = v[:, np.newaxis]
    # print(R.shape)
    # print(v.shape)
    # print(np.matmul(R, v))
    #
    # l = [1, 2, 3]
    # for i in range(5):
    #     try:
    #         print(l[i])
    #     except:
    #         continue

    # image_path = 'C:\\study\\FFR\\myProject\\mytest\\mydata\\input\\prediction_00177.nii.gz'
    # image_path = 'test_rotation.nii.gz'
    # image = sitk.ReadImage(image_path)
    #
    # slice = sitk.GetArrayFromImage(image)
    # slice = slice.swapaxes(0, 2)
    # slice = slice[276]
    # slice = slice.astype(np.uint8)
    # slice = slice[:, :, np.newaxis]
    # contours, _ = cv2.findContours(slice, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_NONE)
    # rect = cv2.minAreaRect(contours[0])
    # print(rect)
    #
    # center = (276, rect[0][0], rect[0][1])
    # trans = sitk.Euler3DTransform()
    # trans.SetCenter(center)
    # trans.SetRotation(np.deg2rad(90 - rect[-1]), 0, 0)
    # resampled = sitk.Resample(image, trans)
    #
    # sitk.WriteImage(resampled, 'test_rotation_x.nii.gz')
    # plt.imshow(sitk.GetArrayFromImage(image)[:, :, 276])
    # plt.show()
    # plt.imshow(sitk.GetArrayFromImage(resampled)[:, :, 276])
    # plt.show()

    # array = np.asarray([[1.0, 0],
    #                     [0, 1]])
    # a = np.logical_or(array, np.asarray([[0, 1.0], [0, 0]]))
    # print(a.astype(np.float32))

    # test rotate to axis
    # np.set_printoptions(precision=16)
    # img_path = r'C:\study\FFR\myProject\mytest\mydata\output_with_axis\outputAxis_00186.nii.gz'
    # axis_path = r'C:\study\FFR\myProject\mytest\mydata\generated_axis\generated_00186.txt'
    #
    # image = sitk.ReadImage(img_path)
    # vec = np.zeros(3)
    # with open(axis_path) as f:
    #     text = f.readlines()
    #     vec_line = text[1].replace('[', '').replace(']', '').split()
    #     vec[0] = float(vec_line[0])
    #     vec[1] = float(vec_line[1])
    #     vec[2] = float(vec_line[2])
    # theta_y, theta_z = geoutils.get_rotation_angle(vec)
    # print(theta_y, theta_z)
    # trans = sitk.Euler3DTransform()
    # # trans.SetCenter((
    # #     image.GetSize()[0] / 2,
    # #     image.GetSize()[1] / 2,
    # #     image.GetSize()[2] / 2
    # # ))
    # # trans.SetRotation(0.0, 0.0, theta_z)
    # trans.SetTranslation((200, 0, 0))
    # out = sitk.Resample(image, trans)
    # # trans.SetRotation(0.0, theta_y - 0.2269, 0.0)
    # sitk.WriteImage(out, 'test_rotate_z.nii.gz')
    # img = sitk.ReadImage('test_rotate_z.nii.gz')
    # trans.SetRotation(0.0, theta_y, 0.0)
    # # trans.SetTranslation((0.0, 0.0, 0.0))
    # out2 = sitk.Resample(img, trans)
    # sitk.WriteImage(out2, 'test_rotate_y.nii.gz')

    # test_vec = np.asarray([1.0, 1.0, 1.0])
    # test_vec /= np.linalg.norm(test_vec)
    # array = np.zeros((320, 512, 512))
    # temp_point = np.asarray(array.shape) / 2
    # start = temp_point.copy()
    # step = 2
    # end = np.asarray([0.0, 0.0, 0.0])
    # for i in range(200):
    #     ind_i = round(temp_point[0])
    #     ind_j = round(temp_point[1])
    #     ind_k = round(temp_point[2])
    #     for j in range(-2, 3):
    #         for k in range(-2, 3):
    #             try:
    #                 array[ind_i, ind_j + j, ind_k + k] = 2.0
    #                 end = np.asarray([ind_i, ind_j, ind_k])
    #             except:
    #                 continue
    #     temp_point += step * test_vec
    # print(start, end)
    # draw_vec = end - start
    # draw_vec /= np.linalg.norm(draw_vec)
    # print('theta between draw_vec and test_vec: ', np.rad2deg(math.acos(draw_vec.dot(test_vec))))
    #
    # sitk.WriteImage(sitk.GetImageFromArray(array), 'test_rotate_correct_origin.nii.gz')
    #
    # out_image = sitk.GetImageFromArray(array)
    #
    # out_image.CopyInformation(image)
    # theta_y, theta_z = geoutils.get_rotation_angle(test_vec)
    # print(theta_y, theta_z)
    # trans = sitk.Euler3DTransform()
    # trans.SetCenter((
    #     out_image.GetSize()[0] / 2,
    #     out_image.GetSize()[1] / 2,
    #     out_image.GetSize()[2] / 2
    # ))
    # trans.SetRotation(0.0, 0.0, theta_z)
    # out = sitk.Resample(out_image, trans)
    # sitk.WriteImage(out, 'test_rotate_correct_z.nii.gz')
    # trans.SetRotation(0.0, theta_y, 0.0)   # have an error of 10 degree, aka 0.1745
    # out = sitk.Resample(out, trans)
    # sitk.WriteImage(out, 'test_rotate_correct_y.nii.gz')

    # np.set_printoptions(precision=16)
    # with open('test_precision.txt', 'w') as f:
    #     point = np.asarray([1.2345678987654321111, 1.2345678987654321, 1.2345678987654321])
    #     f.write(str(point[0]) + ', ' + str(point[1]) + ', ' + str(point[2]) + '\n')
    #     f.write(str(point) + '\n')
    #
    # print(np.get_printoptions())

    # for i in range(5):
    #     try:
    #         a = 1 / None
    #     except:
    #         a = i
    #     print(i)

    # image_path = r'C:\study\FFR\myProject\mytest\test_skele_from_anti.nii.gz'
    # slice = 80
    # image = sitk.ReadImage(image_path)
    # slice_array = sitk.GetArrayFromImage(image)[slice]
    # fig = plt.figure()
    # ax = fig.add_subplot(1, 2, 1)
    # ax.imshow(slice_array)
    #
    # k = morphology.disk(10)
    # res = morphology.closing(slice_array, k)
    #
    # ax = fig.add_subplot(1, 2, 2)
    # ax.imshow(res)
    # plt.show()

    # k_3d = morphology.ball(8)
    # array = sitk.GetArrayFromImage(image)
    # t = time.process_time()
    # out = morphology.closing(array, k_3d)
    # tt = time.process_time()
    # print('==========finish, cost', tt - t, 'seconds=============')
    # out_img = sitk.GetImageFromArray(out)

    # array = sitk.GetArrayFromImage(image)
    # k = morphology.disk(15)
    # t = time.process_time()
    # for z0 in range(array.shape[0]):
    #     array[z0] = morphology.closing(array[z0], k)
    # tt = time.process_time()
    # print('----finish z axis, cost', tt - t, 'seconds----')
    # t = tt
    # for y0 in range(array.shape[1]):
    #     array[:, y0, :] = morphology.closing(array[:, y0, :], k)
    # tt = time.process_time()
    # print('----finish y axis, cost', tt - t, 'seconds----')
    # t = tt
    # for x0 in range(array.shape[2]):
    #     array[:, :, x0] = morphology.closing(array[:, :, x0], k)
    # tt = time.process_time()
    # print('----finish x axis, cost', tt - t, 'seconds----')
    # out_img = sitk.GetImageFromArray(array)
    #
    # sitk.WriteImage(out_img, 'test_3d_x_closing.nii.gz')

    # TODO: test closing before blur
    # image_path = r'C:\study\FFR\myProject\mytest\mydata_apply\input\case_00218\segmentation.nii.gz'
    # image = sitk.ReadImage(image_path)
    # array = sitk.GetArrayFromImage(image)
    # array = array.astype(np.uint8)
    # slice = 339
    #
    # fig, axes = plt.subplots(2, 2, sharex=True, sharey=True)
    # ax = axes.ravel()
    #
    # ax[0].imshow(array[:, :, slice])
    # anti_array = array.copy()
    # ax[1].imshow(cv2.medianBlur(array[:, :, slice], 9))
    #
    # k = morphology.disk(10)
    # closed = morphology.closing(array[:, :, slice], k)
    # k_3d = morphology.ball(10)
    # # closed_3d = morphology.closing(array, k_3d)
    #
    # ax[2].imshow(closed)
    # ax[3].imshow(cv2.medianBlur(closed, 9))
    # # ax[2].imshow(closed_3d[:, :, slice])
    # # ax[3].imshow(cv2.medianBlur(closed_3d[:, :, slice], 9))
    #
    # plt.show()

    # d0 = np.asarray([-0.7921039184295339, 0.6103862567330447])
    # p0 = np.asarray([206., 240.76923076923077])
    # d1 = -1 * np.asarray([0.9682004730594804, 0.25017562624564])
    # p1 = np.asarray([206., 356.0769230769231])
    # pi = np.asarray([79.76452937750912, 331.8141300966583])
    #
    # dd0 = pi - p0
    # dd0 /= np.linalg.norm(dd0)
    #
    # dd1 = pi - p1
    # dd1 /= np.linalg.norm(dd1)
    #
    # print(np.rad2deg(math.acos(dd0.dot(d0))))
    # print(np.rad2deg(math.acos(dd1.dot(d1))))

    # TODO: test rodrigues

    image = sitk.ReadImage('C:/study/FFR/myProject/mytest/mydata_test/images/case_00177/segmentation.nii.gz')
    t = sitk.Euler3DTransform()
    center = np.asarray(image.GetSize()) // 2 * np.asarray(image.GetSpacing())
    t.SetCenter(center)
    R = geoutils.get_rodrigues_rotation_matrix([1., 0., 0.], [0., 0., 1.]).flatten().tolist()
    print(R)
    t.SetMatrix(R)
    out = sitk.Resample(image, t)
    sitk.WriteImage(out, 'test_rodrigues.nii.gz')

    pass
